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Abstract. An axisymmetric collapse of non-rotating gravitational waves is 
numerically investigated in the subcritical regime where no black holes form but 
where curvature attains a maximum and decreases, following the dispersion of the 
initial wave packet. We focus on a curvature invariant with dimensions of length, 
and find that near the threshold for black hole formation it reaches a maximum 
along concentric rings of finite radius around the axis. In this regime the maximal 
value of the invariant exhibits a power-law scaling with the approximate exponent 
0.38, as a function of a parametric distance from the threshold. In addition, 
the variation of the curvature in the critical limit is accompanied by increasing 
amount of echos, with nearly equal temporal and spatial periods. The scaling 
and the echoing patterns, and the corresponding constants are independent of the 
initial data and coordinate choices. 



PACS numbers: 04.25.D-,04.25.dc,04.20.-q 

1. Introduction 

Universality, scaling and self-similarity found in critical gravitational collapse is 
one of the most fascinating phenomena associated with gravitational interactions. 
First discovered numerically by Choptuik [T] in spherically-symmetric collapse of 
massless scalar field, this distinctive behaviour was later observed in other systems, 
including those with various matter contents and equations of state, diverse spacetime 
dimensions etc. However, while a great deal of literature has emerged on critical 
phenomena in spherical symmetry, only a limited number of non-perturbative studies 
exists in less symmetric settings, see [5] for a review. 

Perhaps the simplest non-spherical system is a pure axisymmetric gravitational 
wave, collapsing under its own gravity. Abrahams and Evans [3] found that the mass of 
black holes, forming in the evolution of sufficiently strong initial waves, exhibits scaling 
of the form Mbh oc (a — a») /3 in the limit when the strength parameter a tends to a*, 
the threshold for black hole formation, and determined the exponent of the power-law 
to be (3 ~ 0.35 — 0.38. They have also given less conclusive evidence of periodic echoing 
of the near critical solutions. Surprisingly, these results proved difficult to reproduce; 
in fact no other successful simulation of the axisymmetric vacuum collapse has been 
reported to date (see e.g. [4] for a failed attempt). However, in this paper we present 
new results obtained with the aid of our recent harmonic code [5]. 

We focus on subcritical collapse of axisymmetric non-rotating Brill waves during 
which black holes do not form, but where curvature grows to reach a maximum and 
subsequently diminishes, following the dispersion of the initial wave. Perturbative 
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studies of the critical solutions [5J[S] suggest that the power-law scaling of characteristic 
quantities near the critical point should occur on both sides of the black hole formation 
threshold, regardless of appearance of horizons. While this was confirmed in numerical 
experiments in spherical symmetry, see e.g [71 [5] , it is an open question whether same 
is true in other situations as well. Here we demonstrate that in the axisymmetric 
subcritical collapse, a curvature invariant with dimensions of length follows a power- 
law with the exponent, j3 = 0.385 ± 0.015, similar to that found by Abrahams 
and Evans in supercritical case. Additionally, we find that the solutions develop 
increasingly large numbers of echos as the critical limit is approached. Our current 
resolution allows observation of up to three echos around the time instant where 
curvature is maximal. We measure that, for example, the Riemann curvature invariant 
oscillates in time with the (logarithmic) period of A T ~ 1.1, and that the logarithm 
of the invariant changes on each echo by nearly same amount A r ~ A T ~ 1.1. 

We verify that the scaling and echoing constants are essentially independent of 
initial data and specific coordinate conditions used to calculate the solutions. In 
contrast to spherically symmetric collapse, where the greatest curvature is always 
at the origin, the evolution of the axisymmetric waves is more complicated and the 
spacetime location of the maximum depends strongly on the geometry of the initial 
data. We evolved series of subcritical initial data where curvature attained a maximum 
along equatorial rings of various radii centered at the axis. Besides, we found that 
in supercritical evolutions of the same data an apparent horizon forms, engulfing the 
ring-shaped locus of the maximal curvature. This indicates that the critical solutions 
found in the Brill-wave evolutions are different from ones calculated by Abrahams & 
Evans, in whose case the maximal curvature has always occurred at the origin, and 
the black holes tend to be arbitrary small in the critical limit. Strikingly, despite these 
differences the near critical scaling and echoing patterns are similar, and the scaling 
exponents are comparable. 

In the next section we briefly describe the initial value problem for constructing 
the axisymmetric vacuum asymptotically flat spacetimes without angular momentum; 
the details of the equations, gauge conditions and our numerical code are found in [5]. 
Section [3] is devoted to the results and numerical tests. We summarize our findings, 
discuss limitations of the current method and outline perspectives in the concluding 
section |U 

2. The equations and a method of their solution 

We are interested in solving the vacuum Einstein equations 



where is the Ricci tensor. We consider axisymmetric asymptotically flat 

spacetimes without angular momentum, and assume they can be foliated by a family 
of spacelike hypersurfaces, starting with the initial surface at t = 0, where the spatial 
metric and its normal derivatives are chosen to satisfy the constraints (Gauss-Codazzi 
equations). 

The most general metric adapted to the symmetries of the problem can be written 
using the cylindrical coordinates 



Rni> = o, 



(i) 



ds 2 = g ab dx a dx b + r 2 e 2S d(f> 2 , 



(2) 
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where the seven metric functions — g a j,, a, b = 0, 1, 2 and S — depend only on t, r and z. 
[3 In order to solve the field equations ((T|) we employ the Generalized Harmonic (GH) 
formalism [SJ [TUJ QT] , adopted to the axial symmetry in [5] . To this end we define the 
GH constraint, 

C a = -Dx a + H a = ~T a afj g af} + H a = 0, (3) 

where r°g are Christoffel symbols, and the "source functions" H a = H a (x,g) depend 
on coordinates and the metric (but do not on the metrics's derivatives) and are 
arbitrary otherwise. We than modify the Einstein equations: 

Rye, - C M = 0, (4) 

that now become a set of quasilinear wave equations for the metric components of the 
form g a P g[_iv,at3 + . . . = 0, where ellipses designate terms that may contain the metric, 
the source functions and their derivatives. 

Fixing the coordinate freedom in the GH language amounts to specifying the 
source functions, and we choose those by requiring that the spatial coordinates satisfy 
damped wave equations, while the time coordinate remains well behaved while the 
lapse satisfies a damped wave equation |12|I5]. A particular example of these conditions 
[12"] . that we use here, can be written in terms of the kinematic ADM variables as 

H° w = 2fi x \og(^-\ n a -2[i 2 a- 1 lal (3\ (5) 



where — (—g 00 )~ 1 ^ 2 d IJ- 1 is the unit normal to the spatial hypersurfaces of constant 
time, a is the lapse, f3 l is the shift, y a b — g a t + n a nt is the spatial metric, 
7 = (511 1722 — 512) ex P(S), and fix and [12 are parameters. 

The initial data is given at t = 0, where we choose the initial spatial metric to be 
in the form of the Brill-wave [13] 

dsl = V> 4 (r, z) \ e 2rB ^ (dr 2 + dz 2 ) + r 2 dcjA , (6) 

with 

( r 2 z 2 \ 

B(r,z)=arexp[ , (7) 



where ay, o~ z and a are parameters. 

We further assume time-symmetry, in which case the momentum constraint 
identically vanishes at t = 0, while the Hamiltonian constraint becomes the elliptic 
equation for ip 

d 2 + l-d r + d^j i, = -^r (d 2 r + \ + d%j B, (8) 

which is solved subject to regularity conditions at the axis, equatorial reflection 
symmetry, and asymptotic flatness boundary conditions: 

d z ip{r,0)=0, d r ip(0,z) = 0, ip(r, 00) = ^(oo,z) = 1. (9) 
We assume initially harmonic coordinates, H a = 0, and choose the initial lapse 
a(t = 0, r, z) = gl(?(t = 0, r, z) = 1. 

Having specified the initial data we integrate the equations Q forward in time, 
imposing asymptotic flatness and regularity at the axis, r = 0. For simplicity, 
we restrict attention to the spacetimes having equatorial reflection symmetry. The 
highlights of our finite-differencing approximation (FDA) numerical code [5] that we 
employ to solve the equations include: 

I While Greek indices range over t, r,Z,(f> — 0, 1, 2, 3, Latin indices run over 0, 1, 2. 
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• An introduction of a new variable that facilitates axis regularization. While 
elementary flatness at the axis implies that each metric component has either 
to vanish or to have vanishing normal derivative at that axis, requiring absence 
of a conical singularity at r = results in the additional condition: gn(t, 0, z) = 
exp[2 S(t, 0, z)]. Therefore, at r = we essentially have three conditions on the 
two fields S and g\\. While in the continuum, and given regular initial data, the 
evolution equations will preserve regularity, in a FDA numerical code this will be 
true only up to discretization errors. Our experience shows that the number of 
boundary conditions should be equal to the number of evolved variables in order 
to avoid regularity problems and divergences of a numerical implementation. We 
deal with this regularity issue by defining a new variable 



that behaves as A ~ 0(r) at the axis, and use it in the evolution equations 
instead of S. This eliminates the overconstraining and completely regularizes the 
equations. Crucially, the hyperbolicity of the GH system is not affected by the 
change of variables. 

• Constraint damping. The constraint^ equations, C M = 0, are not solved in the 
free evolution schemes like ours, except at the initial hypersurface. While one 
can show that in the continuum the constraints are satisfied at all times, in FDA 
codes small initial violations tend to grow and destroy convergence. A method 
that we use to damp constraints violations consists of adding to the equations Q 
the term of the form [HI [16] , 



where k is a parameter. We note that Z^ v contains only first derivatives of the 
metric and hence does not affect the principal (hyperbolic) part of the equations. 

• A spatial compactification is introduced in both spatial directions by transforming 
to the new coordinates x = x/(l + x), x G [0,1], x G [0,oo), where x stands for 
either r or z. The advantage of this scheme is that asymptotic flatness conditions 
9nv = vffu a t the spatial infinity are exact. 

• We use Kreiss-Oliger-type dissipation in order to remove high frequency 
discretization noise. [| An additional role of the dissipation is to effectively 
attenuate the unphysical back reflections from the outer boundaries, resulting 
from the loss of numerical resolution there. This allows using compactification 
meaningfully [TTJ [17] . 

In order to characterize the spacetimes that we construct, we use the Brill mass 
|13j . computed at the initial time-slice, 



which — we verify — coincides with the ADM mass. For the purpose of quantifying the 
strength of the gravitational field we calculate the Riemann curvature invariant having 
dimension of inverse length, 




(11) 




(12) 



i = (R a ^R a ^n 1/4 



(13) 



§ It can be shown that the standart Hamiltonian and momentum constraints are equivalent to the 
GH constraints |14| . 

| i.e. the noise with frequency of order of the inverse of the mesh-spacing 
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0.2, 1.1 
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0.04 


8.21 ±0.01 


0.593 


0.0 
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Table 1. The parameters of the initial data ov,cr z , as well as grid spacing 
h and gauge parameters pi,p.2 determine the threshold amplitude a* whose 
upper margin corresponds to formation of a black hole, and whose lower margin 
corresponds to a regular spacetime. Given this set of initial parameters, this 
further determines the total mass M* and the "accumulation locus", whose 
position and time of occurrence is given by the radial position p* and proper 
time t, . The radial position p„ is measured in terms of the circumferential radius 
1151 . and the proper time t* at that location is measured in units of the total 
mass. The parameters 21 and 41 indicate that 2 and 4 AMR levels were used, 
respectively. All other simulations are unigrid. 



at various locations, and in certain experiments we follow its evolution in the proper 
time at that location (r, z), 

r(t,r,z) = [ a{t' ,r,z)dt' . (14) 



ii 



We also use the circumferential radius 

p = re s . (15) 



3. Results 



The initial data (16I7[) are characterized by the amplitude a and the "shape" parameters 
ay and a z , which define the mass of the data and their "strength" , namely the tendency 
to collapse and form a black hole. For a given amplitude and fixed oy + a z = const, 
the data with larger a z /a r are stronger (see also [T^l HH|)- In addition, by varying 
the shape parameters at fixed gauge, we can control the spacetime locations where 
curvature evolves to a maximum or where an apparent horizon first forms. In our 
experiments we use several sets of ay and o~ z , and adjust the strength of the initial 
wave by tuning its amplitude. 

The initial data are numerically evolved forward in time. We use grids with similar 
mesh-sizes in both spatial dimensions h r = h z = h, and time-steps of h t = 0.04 h and 
h t = 0.05 h. Usually our fixed grids consist of 200, 250, 300 or 400 points, uniformly 
covering the compactificd spatial directions. We also experiment with adaptive mesh 
refinement (AMR), provided by the PAMR/AMRD software [2D] . In this case we use 
two or four refinement levels, and the base mesh with the resolution of ft. = 1/128. The 
Kreiss-Oliger dissipation parameter is typically cko ~ 0.5 — 0.85, with larger values 
used on finer grids and stronger initial data; and the constraint damping parameter 
in (|TT|) is k — 1.4 — 1.7. The gauge fixing parameters ([5]) in the ranges /ii ~ 0.1 — 0.3 
and /x 2 ~ 0.9 — 1.2, usually gave stable, sufficiently long evolutions. 



On critical collapse of gravitational waves 



6 



The system is weakly gravitating for small amplitudes, in which cases the initial 
wave packet ultimately disperses to infinity. However, for amplitudes above certain 
threshold, a*, the wave collapses to form a black hole, signaled by an apparent horizon. 
In subcritical spacetimes we can define the "accumulation locus" where curvature 
attains a global maximum before decaying. In our coordinates ([5]), and for our initial 
data (where the ratio of cr's never exceeds five) the position of the maxima (t*,r*, z*) 
is always along the equator z* = 0. 

The threshold amplitude for black hole formation, a*, depends on the initial 
data, controlled by a r ,a z , and gauge parameters /ii,/i2, and the resolution, h. Table 
Q] records critical amplitudes, a*, masses, M* and the spacetime positions, /9*,t», of 
the accumulation locus in the strongest, a ~ a*, initial data evolutions, for a few 
sets which we have calculated. In contrast to spherically-symmetric collapse, where 
accumulation locus is solely at the origin, in axial symmetry this is not always the case. 
For instance, the critical amplitude for the initial data with a r = a z = 1, determined 
in the unigrid simulations with h = 1/300 is a* = 6.20021. The spacetime position of 
the accumulation depend on the amplitude such that for a J$ 0.99a* the accumulation 
loci are at the origin, and for larger amplitudes they shift to be along the rings of radii 
p, ~ 0.2. The time of occurrence of the maxima converges to t» ~ 1.46 M* from above 
in the limit a — > a* . While qualitatively similar behaviour is observed for most of the 
initial data families listed in TableQ] the initial data defined by ay = 0.7, a z = 1.5 has 
the accumulation loci at the origin all the way to the strongest subcritical amplitude 
of a = 8.20. However, since in this case we have succeeded to compute a» only with a 
modest accuracy of 1 part in 820, a possibility remains that closer to the threshold the 
accumulation loci will become ring-shaped. For this set the time of the accumulations 
converges to r* ~ 3.2 M* in the limit a — > a*, and the mass of the near critical 
solutions, M* ~ 0.593, is about one half of that found in the ay = 0.9, a z = 1.3 case. 

As described next, there is a power-law scaling of the maximal curvature in the 
limit a/a* — > 1 for all families of the initial data listed in Table [TJ In most cases 
the scaling shows up at relatively large values of a*/a — 1 ~ 10~ 3 , for all resolutions 
better than h = 1/200. However, it turns out that the data calculated in fixed-mesh 
simulations with h > 1/250 is too noisy and dependent on the details of numerics, to 
provide a reliable estimate of the scaling exponent. 

The scaling can be envisaged by plotting the maximal value of the Riemann 
curvature invariant (|13[) as a function of the parametric distance from the critical 
amplitude, a* — a; this is shown in Fig. [TJ Each point here represents the global 
maximum |/ max | computed during evolutions defined by ay = 0.9 and o~ z = 1.3, and 
the numerical parameters: h = 1/300, h t /h = 0.04, fix = 0.3, = 0.9, k = 1.7, e = 0.6. 
The solid line represents the least-squares linear fit to the data. The slope of the line, 
/3 ~ —0.37, is in agreement with the exponent of the black- holes' mass scaling, found 
in supercritical collapse by Abrahams and Evans [3]. The data depicted in Fig. [2] were 
obtained with a differently shaped initial wave, a r — a z ~ 1, and the parameters 
h = 1/300, ht/h = 0.05, Mi = 0.12, ^ 2 = 1.17, K = 1.7, e = 0.8. The threshold 
amplitude in this case is found with somewhat lesser accuracy, a» = 6.20021 ±0.00001. 
However, the data is still fitted well with a straight line whose slope, /3 ~ —0.4, 
coincides with the exponent in Fig. [ljto within 8%. 

K For comparison, in scalar field collapse the signatures of near-critical scaling do not appear before 
a,/a - 1 < 10" 8 . 

+ Note that our exponent is negative since the dimensions of / are inverse length, while black hole 
mass computed in [3] has dimensions of length. 
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Figure 1. A logarithm of the maximal Riemann invariant H13I I as a function 
of the distance from the critical amplitude, a, — a, in the simulations with 
o r = 0.9, <J z = 1-3, and the fixed resolution h = 1/300. The critical amplitude 
in this case is a„ = 7.3079082, and the maximal curvature is \I m ax\ ~ 10 5 in 
the units of the total mass. The linear fit to the data (solid line) has the slope 
~ —0.37. Notice the (quasi-) periodic "wiggle" of the data points about the 
straight line, which we interpret to signal periodic self-similarity of the critical 
solution. 
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Figure 2. A plot similar to Fig. [l] but obtained with different parameters: 
cr r = <Tz = 1 and the resolution h = 1/300. In this case the critical amplitude is 
a, = 6.20021. Remarkably, the slopes of the linear fits in both figures agree to 
within 8%. 

It is remarkable that despite the evolutions of the initial waves shown in Figs. 
[T] and [2] are dramatically different, the maximal curvatures in both cases follow a 
power-law with similar exponents. We verify that the same scaling again appears in 
simulations with other shape parameters and in all cases that resulting exponent is 
consistently in the range j3 ~ 0.37 — 0.4. In addition, the scaling exponent within these 
bounds when we use coordinate conditions with different choices of /i's in ((SJ) (see e.g. 
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Figure 3. The dynamics of the logarithm of the Riemann invariant, /, as a 
function of the proper time r(r») for several values of a/a*, obtained in the 
evolution of initial data defined by u r = 0.9, <r z = 1.3, and the fixed resolution 
h = 1/300. The variation of I toward the accumulation, r„, and away from it is 
accompanied by oscillations, whose number grows in the limit a — > a„ . The double 
dip in top right panel at about t ~ 1.51 and in bottom panels near r ~ 1.48 is 
a result of the interference between the main and a secondary reflection off the 
axis. 



Figs. [I] and [2]). While this does not test the rigidity of (3 with respect to all possible 
coordinate conditions, this demonstrates relative consistency of the exponent within 
the large family of the gauges ([5]) . We conclude that in the critical limit the maximal 
curvature predominantly scales as |/ m aa;| oc (a* — a) - ' 9 , with /? = 0.385 ± 0.015, where 
the errorbars represent the deviation from the average value computed over all initial 
data sets that we have evolved0 

The distribution of data points in Figs. [1] and [2] has a striking property, namely 
the data "wiggle" about the linear fit. We note that similar wiggle was also observed 
in near critical collapse of scalar field. In that case it was attributed [52] to the 
periodic self-similarity found in that system, where the critical solution, Z*, repeats 
on itself after a discrete period A: Z*(t, r) — Z*(t e A ,r e A ). Besides, |H] found 
that the period of the wiggle is A/(2 /3), and thus may, in principle, allow calculating 
the self-similarity scale A by measuring the slope and the period of the wiggle in a 
plot like ours Figs. Q] and [2] We believe that the quasi-periodic fluctuations of the 
points about the linear fit in these figures do signal discrete self-similarity, however 
our current data are insufficiently accurate and have too short a span to provide more 
quantitative estimate of the wiggle period, beyond a very rough value of anything 
between two and four. 

* The slope obtained for each initial data set carries individual fitting errors. However, these are 
typically smaller than the fluctuations around the average /3 computed over all data sets. 
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Figure 4. The echoing pattern obtained in the evolution of the initial data sets 
with o r = 0.9, cr z = 1.3 and a = 6.940. Left panels show the low resolution runs 
that use 2 levels of AMR, with the base h = 1/128, other panels were obtained 
using 4 AMR levels with the same base; the rightmost panels is the zooming of 
the late time behavior shown in middle panels. While the lower resolution runs 
diverge around r ~ 1.47, the higher resolution runs extend beyond that, allowing 
to calculate additional cchos. Notice that / has sharper and easier to identify 
features than a. 

Independent, and more direct signatures of discrete self-similarity are obtained 
by examining the behaviour of the curvature when a — > a„. It turns out that in this 
limit, in addition to that / attains increasingly larger maxima, the temporal variation 
of / is also accompanied by increasing amount of oscillations. This is illustrated in 
Fig. [3l which shows the variation of / as a function of the proper time, calculated 
at the accumulation loci, for a sequence of a's. Figure shows that the amount of 
fluctuations — indicated by the peaks or inflection points — grows from one to three 
in the limit a/a* — > 1 on both sides of the accumulation locus. Such an oscillatory 
behaviour is again reminiscent of the "echoing" in critical spherical collapse of scalar 
field (see e.g. Fig. 5 in [8] and Fig. 7 in [21] ) and we interpret it as evidence of 
periodic self-similarity in our system as well. 

Like the power-law scaling of the maximal curvature, the echoing of our solutions 
in the near critical limit is independent of specific gauges or particular initial data sets. 
This is demonstrated in Fig. |4l which depicts temporal evolutions of the Riemann 
invariant and the lapse function found in simulations of the initial data characterized 
by ay = 0.9, a z = 1.3, the gauge constants /ii = 0.3, [i2 = 0.9, and the amplitude 
a = 6.940. The functions in left panels were computed using 2 levels of AMR, and 
the other panels were generated using 4 levels of AMR; in both cases the base-level 
resolution is h = 1/128. Figure shows that the dynamics in this case involves more 
scatterings and interferences of the initial and secondary waves than e.g. in ay = ay 
runs, depicted in Fig. [3] 

In most cases higher resolution simulations run longer and allow computing more 
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Figure 5. The typical temporal variation of the curvature invariant I near the 
accumulation locus is oscillatory in time. Shown is the evolution of the initial 
data with a r = <r z = 1 and a /a* = 0.999998. On each oscillation log |/| varies by 
A r ~ 1.1 ± 0.1, which is close to the time period A T ~ 0.95 ± 0.15 of the four 
oscillations around t*. 



oscillations. Notice that the shapes of the curves in left and middle panels in Fig. 
[4] are essentially identical until r ~ 1.47. However, while the lower resolution runs 
diverge around that time due to formation of a singularity, the higher resolution 
evolutions continue beyond that, and develop additional echos that accumulate near 
t* ~ 1.477M*, just before the numerics fails. The critical amplitude, determined in the 
4-level AMR simulations is a» = 6.9401, and the accumulation loci occur are the origin. 
While we were unable to stabilize the 4-level evolutions for amplitudes beyond about 
a ~ 7, in lower-resolution, 2-level runs we find a different critical solution with the 
amplitude, a* = 7.246067, where the accumulation locus lies at p* ~ 0.23, see TableQ] 
The total masses of the near critical spacetimes, the accumulation loci and such details 
of evolutions as the amount of secondary scatterings and interference's, reflected in 
the strong variability of the curvature profile, and the total amount of gravitational 
radiation are different in near-critical evolutions in the 2 and 4-level AMR simulations. 
Nevertheless, the scaling and echoing constants appear to be nearly identical. 

In order to estimate the period of the echos we plot in Fig. [5] the temporal 
variation of I computed in the evolution of the initial data set having er r — a z = 1 
and a = 6.2002. By measuring the distances between the peaks or inflection points — 
marked by arrows in Fig. [5] — we find that the curvature fluctuates in time with the 
logarithmic period A r = 0.95 ± 0.15 and that on each echo the logarithm of / varies 
by approximately A r ~ 1.1 ±0.1. The errorbars here represent the maximal deviation 
from the average values of A T and A r , measured in this figure. We note that both 
periods agree within the error-bars. A similar Fig. [6] shows the dynamics of / against 
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Figure 6. The variation of I obtained in a 4-level AMR simulations initialized 
by oy = 0.9, cr z = 1.3, and a/a* ~ 0.99997. The variation of log|/| on each 
oscillation (marked by arrows) is nearly equal to the temporal period of the echos, 
A r ~ A T ~ 1.1. 



r* — t, that was obtained in simulations with 4 levels of AMR, ay = 0.9, a z = 1.3, and 
a — 6.940, shown in right panels in Fig. 2) Although the resulting dynamics is quite 
complicated, featuring multiple scatterings and interferences, there are 3 prominent 
peaks — marked by the arrows in Fig. [B] — that can be identified as echos. Their 
temporal period is A T = 1.10 ± 0.04, and on each echo the logarithm of |/| grows 
by a comparable amount A r = 1.12 ± 0.06. We note that these values match within 
the errorbars, and are in a good agreement with the periods computed in Fig. [5] 

The echoing is not specific to the curvature invariant /, other metric functions 
oscillate as well. However, while the echoes of / are signalled by the sharp peaks 
the fluctuations of metric components are typically milder, showing up as inflection 
points (see bottom panels in Fig. [4]). This makes / a superior quantity for the 
purpose of measuring the echoing periods. Although we mainly discussed variations 
of / at the location of its global maximum, we verified that curvature develops 
echoes in other locations as well, where, however, the amount of the echos and their 
amplitude is generally smaller than around the accumulation locus. Since away from 
the accumulation locus the curvature remains bounded in the critical limit we expect 
only a finite number of such oscillations. 

We conclude this section by briefly discussing the accuracy of our code. While it 
was not possible to perform a direct convergence test of e.g. the critical amplitude or 
the scaling exponent, since changes of the resolution usually required readjustments of 
the dissipation, cko, and the constraint damping, k and the gauge parameters ^1,^2, 
which alter the "conditions of the convergence study" , that require all parameters 
except h to stay fixed. Nevertheless, as indicated in Table Q] the critical amplitude 
seems to converge as a function of the resolution, at least in the equal-cr's case. In 
addition, formal numerical convergence tests in individual, fixed amplitude simulations 

J The evolution of this initial data diverges soon after the accumulation at about r ~ 1.477M*, due 
to imperfections of our AMR numerics, and sensitivity to the choices of fii and 112 ■ Hence only the 
collapse stage of the evolution is shown in this figure. 
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Figure 7. The logarithm of \I m ax\ obtained in simulations with two distinct 
sets of the damping and dissipation constants: k = 1.4, tKO = 0.75 (right panel), 
and k = 1.5, ejfo = 0.85 (left panel). We use same resolution of h = 1/300, and 
all other equal parameters. The difference in /3's in this case is less than 3%, 
indicating the quality of our numerics. 

along with the demonstration that the Hamiltonian and the momentum constraints 
are satisfied during the evolutions were carried out in [5], indicating nearly second 
order convergence, and exponential decay of the ^-norms of the constraints at late 
times. 

The consistency of the scaling exponents obtained in simulations with a whole 
different set of parameters (see e.g. Figs. [Hand [2]) indicates robustness of /?, however 
the overall accuracy of our code can be estimated by changing only the numerical 
parameters. To this end we performed simulations with different damping and the 
Kreiss-Oliger dissipation constants k and zko-, and otherwise similar parameters. Fig. 
[7] shows that /3's computed in two sets of simulations defined by K — 1.4, exo — 0.75, 
and k = 1.5, cko = 0.85 differ by less than 3%. 

4. Discussion 

The ring-shaped accumulation loci that we observe in most evolutions of the time- 
symmetric Brill- wave initial data indicate that the critical solutions in these cases are 
genuinely different from those found by Abrahams & Evans [3] for the ingoing / = 2 
quasi-linear wave data, where the maximal curvature occurs at the origin. The black 
hole radii and their masses found in slightly super-critical evolutions in [3] are tend to 
zero in the critical limit, signalling so-called type II critical phenomenon, characterized 
by smooth transition between dispersion and black hole formation, see Ref. [2]. In 
our case, however, the radii of the accumulation loci are finite, as are the apparent 
horizons that form in supercritical collapse, engulfing the accumulation loci. Although, 
at present, neither our numerics is capable of finding apparent horizons very close to 
the threshold, a^^a^ 1.2 a*, nor it allows to trace the evolution of the horizons to 
their endstate, it indicates that our critical solutions include quasi-stationary ring- 
shaped formations of finite size and mass. 

The situation with the four and two levels AMR simulations is somewhat puzzling. 
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While the two level simulations are clearly divergent near a = 6.9401, which is 
determined as the critical amplitude in the four level runs, comparable in resolution 
unigrid runs do not encounter any particular difficulties at this amplitude. We believe 
such a behaviour may signal another, different critical solution, which is not resolved 
by the lower-dimensional unigrid simulations, and which destabilizes the less accurate 
two-level runs. However, whether this is indeed the case requires further investigation. 

In all cases, we found strong evidence that in subcritical non-rotating 
axisymmetric vacuum collapse curvature exhibits a power-law scaling as a function 
of parametric distance from the threshold for black hole formation. We numerically 
evolved several sets of initial Brill- waves defined by fixed ay and o~ z , and by a tuneable 
amplitude, a, and checked that in the limit a — > a*, |imax| ex (a* — a) ^ with 
roughly the same exponent as that computed in supercritical regime by Abrahams 
and Evans [3J, see Figs. Q]and[2J This demonstrates that quantities with same length 
dimensions — such as the black hole mass in [3] and the inverse curvature invariant 
Imax nere — scale identically. We verified that the exponent is relatively insensitive 
to coordinate conditions. Since we find that the scaling occurs around a ring-shaped 
accumulation locus, which is different from the point-like one of [3], there is no a priori 
reason to expect the exponents in both cases to match. However, the exponents agree, 
and this, apparently, indicates that f3 ~ 0.35 — 0.4 is truly universal and independent 
of the initial data, regardless of what critical solution these data may lead to. [J3 

There is evidence that the near critical solutions are periodically self-similar. 
Specifically, we observe that in the limit a —> a„ the curvature invariant / 
undergoes increasingly large number of oscillations, whose period in the proper 
time is approximately equal to the rate of variation of the curvature on each echo 
A r ~ A r ~ 1.1, see Figs. |3J [5] and [5] We note that the echoing periods reported in 
[3], A ~ 0.6, differ from ours, which are roughly twice larger in magnitude. However, 
this is probably not too surprising since our critical solutions are different from theirs, 
and besides the period of any specific quantity will typically depend on the particular 
combinations of the metric and derivatives, comprising it (for instance, the quantity 
d 2 ^>/dr 2 is twice more variable than "J/). An independent, if circumstantial, signature 
of discrete self-similarity is the distinctive "wiggle" of the data points about the leading 
power-law scaling of \I ma x\, see Figs. [T]and[2] since exactly this kind of behavior is 
expected in the periodically self-similar systems )22j . 




An obvious limitation of the current simulations is their maximal resolution. Even 
though a relatively moderate numerical resolutions of h ~ 1/250 — 1/1000 have already 
provided fruitful insights into the critical behavior, higher resolutions are needed in 
order to compute the scaling and echoing constants more accurately. We expect that 
much closer approach to threshold will be required. This should create a longer span 
of data, enabling a greater accuracy of linear fits in the plots like Fig. Q] and [21 which, 
in turn, will allow unambiguous computation of j3 and of the wiggle period. A closer 
approach a — > a, should also multiply the number of the echoes, allowing a better 
estimate on their periods. Clearly, using numerical meshes of fixed size is not practical 
for probing the limit a — > a*, rather AMR approach should be used. While we have 
already experimented with that, our runs often develop premature instabilities since 
in the near critical limit the system tends to be extremely sensitive to numerical and 

ff In this regard it is interesting to observe that the critical exponent originally found by Choptuik 
in scalar-field collapse, ftsF — 0.374, is again comparable to what we find here. While this may be 
just a coincidence, it may, alternatively, point to the genuine role of the gravity, rather than matter, 
in critical behaviour in scalar-field collapse. 
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gauge parameters; for instance, in the 4-level simulations a slight variation of k by a 
mere 1% ruins convergence. We are currently improving our code in order to locate 
the optimal parameter settings, which will enable us to edge the critical limit, the 
results of that study will be reported elsewhere. 
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